#
# ==============
# Figure A2: Most improved IP rate programs.
# ==============
#
  # All improper payment data.
  DT = fread(file.path(DATALOC,"payments_all.csv"))
  DT[,outlaysamountm := as.numeric(outlaysamountm)]
  DT[,improperpaymentamountm := as.numeric(improperpaymentamountm)]

  # Fix 2021 rates: IP rates for 2021 are in percent, while in pre-2021 in proportions.
  DT[,improperpaymentrate := improperpaymentamountm/outlaysamountm]
  # Drop HUD 2021, which is off.
  DT = DT[!(fiscalyear == 2021 & agency == "hud"),]

  # Limit to programs with at least four observations.
  PT = copy(DT)
  PT[,improperpaymentrate := 100*improperpaymentrate] # convert to percent
  PT[,num_obs := sum(!is.na(improperpaymentrate)),by=c("programid")]
  PT = PT[num_obs >= 4 & !is.na(improperpaymentrate),]
  PT[,program_outlays := sum(outlaysamountm,na.rm=T),by="programid"]

  # Loop over programs and estimate ols slope in time.
  abs = NULL
  for (i in PT[,unique(programid)]) {
    jed = PT[programid == i,lm(improperpaymentrate ~ fiscalyear)]
    abs = rbind(abs,data.table(programid=i,a=jed$coef[1],b=jed$coef[2]),use.names=TRUE)
  }
  # Merge agency and program name for plotting.
  abs = merge(abs,unique(PT[,c("programid","agency","programname","program_outlays")],by="programid",all.x=T))
  setkey(abs,b)
  print(abs)

  # Plots.
  # Medicare Part C.
  f = file.path(FIGURELOC,"FigureA02-PanelA.pdf")
  pdf(f,width=10,height=6)
  par.old = par(mar=c(3.1,3.1,4.1,0.1),cex.lab=1.1,cex.axis=1.1)
  i = 26
  PT[programid == abs[i,programid],plot(x=fiscalyear,y=improperpaymentrate,type="b",pch=19,ann=F,axes=F)]
  abs[i,abline(a,b,lwd=3)]
  x_at = PT[programid == abs[i,programid],seq(min(fiscalyear),max(fiscalyear))]
  axis(1,at=x_at,label=x_at); axis(2,las=2)
  abs[i,title(main=paste(toupper(agency),"\n",toupper(programname),"\nTotal outlays = $",prettyNum(program_outlays,","),"million"))]
  par(par.old)
  dev.off()
  
  # Medicare Part D.
  f = file.path(FIGURELOC,"FigureA02-PanelB.pdf")
  pdf(f,width=10,height=6)
  par.old = par(mar=c(3.1,3.1,4.1,0.1),cex.lab=1.1,cex.axis=1.1)
  i = 31
  PT[programid == abs[i,programid],plot(x=fiscalyear,y=improperpaymentrate,type="b",pch=19,ann=F,axes=F)]
  abs[i,abline(a,b,lwd=3)]
  x_at = PT[programid == abs[i,programid],seq(min(fiscalyear),max(fiscalyear))]
  axis(1,at=x_at,label=x_at); axis(2,las=2)
  abs[i,title(main=paste(toupper(agency),"\n",toupper(programname),"\nTotal outlays = $",prettyNum(program_outlays,","),"million"))]
  par(par.old)
  dev.off()
